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1.0 SUMMARY 


A finite difference method for solving the unsteady transonic flow about harmonically 
oscillating wings is investigated. The procedure is based on separating the velocity potential 
into steady and unsteady parts and linearizing the resulting unsteady differential equation 
for small disturbances. The differential equation for the unsteady velocity potential is linear 
with spatially varying coefficients. 

The work of this report is a direct extension of earlier studies and includes the development 
and application of an out-of-core direct solution. 

The main results of this study are as follows : 

1 . An out-of-core solution procedure has been developed and programmed to provide a 
capability for solving two-dimensional problems with a sufficient number of mesh 
points for practical use. 

2. Reasonable correlation with more exact linear theory is obtained for airfoils of 
vanishing thickness at values of Mach number and reduced frequency of direct 
interest in flutter analyses. 

3. Comparison of two-dimensional finite difference solutions with exact analytic 
solutions indicates the accuracy of the difference solution is dependent on the 
boundary conditions of the mesh region. A corresponding study with similar results 
for a simplified one-dimensional problem is described in NASA CR-2933. 


2.0 INTRODUCTION 


The purpose of the work described in this report is to continue the development of a means 
for calculating air forces for use in flutter analyses of three-dimensional lifting surfaces in 
the transonic flight regime. The work concentrates on a particular procedure that assumes 
small perturbations, the existence of a velocity potential and simple harmonic motion, and 
uses finite difference theory to solve the resulting set of partial differential equations. This 
study represents a direct extension of the research described in references 1 through 3. 

Recent papers by Tijdeman (ref. 4) and Ashley (ref. 5) have contributed significantly 
toward understanding the transonic characteristics that are important to the flutter 
problem. The presence of a shock on an oscillating airfoil results in a pressure pulse in the 
unsteady pressure distribution. This pressure pulse or “pressure doublet” results from the 
chordwise movement of the shock. The shock moves chordwise as the airfoil oscillates. 

The magnitude of the pulse varies as the airfoil moves, and the pulse itself may or may not 
move chordwise. The important characteristics of the pressure pulse with respect to flutter 
instabilities appear to be its amplitude, its phasing with respect to the section motion, and 
its location relative to the section elastic axis. In particular, for airfoils with blunt leading 
edges such as supercritical sections, the pulse tends to be forward on the section, which 
increases the likelihood of instabilities. 

The procedure of this report divides the velocity potential into steady and unsteady parts. 
The steady potential is calculated using the classic, nonlinear, small perturbation differential 
equation. The unsteady potential is then calculated using a linear equation with spatially 
varying coefficients that depend on the steady flow. The procedure does not appear to 
include shock motion since the difference operator that is used for each mesh point is deter- 
mined by the steady flow and does not change during the unsteady calculation. 

The results calculated using this procedure do contain the pressure pulse that results from 
the shock motion (refs. 1 through 3). Figure 1 shows that the pressure coefficient distribu- 
tion computed by subsonic linearized theory for a Mach number of 0.9 indicates no influence 
of a shock. The real and imaginary parts of the pressure jump across an NACA 64A006 air- 
foil are plotted in figure 2 at the same Mach number by the transonic method discussed here. 
By comparing the location of the sharp peaks in figure 2 with the shock in figure 3, we 
readily see that these sharp peaks appear at the location of the steady-state shock on the 
airfoil and hence may well be caused by shock motion induced by the oscillating flow. 

Since we have assumed harmonic motion (i.e., e lcot ), the shock motion could only be repre- 
sented to the first harmonic at best. 

A major concern at the beginning of the investigation was the lack of correlation between 
the finite difference results for a purely subsonic configuration and more exact linear theory. 
Section 6.0 is devoted to this problem, and encouraging results are presented. 

The effect of the choice of boundary conditions on solution accuracy for the two- 
dimensional problem is investigated in section 7. This study parallels an investigation of 
the one-dimensional problem presented in reference 3. 
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Figure 7.- Jump in Pressure Coefficient Across a Flat Plate Oscillating in Harmonic Pitch, 
M = 0.90, co = 0.06 
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Figure 2.— (Concluded) 
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NACA 64A006 airfoil at zero angle of attack 


Data from reference 2 
Air ( 7 = 1.4) 

Freon ( 7 = 1.135) 


Figure 3.— Steady-State Pressure Coefficient for an Airfoil Section; M = 0.9 


Section 8.0 briefly describes several auxiliary investigations that are presented in reference 
6. Finally, conclusions are presented in section 9.0. 

Other results and studies associated with work presented in this document are given in 
references 6, 7, and 8. Auxiliary analytical studies are included in reference 6. Reference 7 
is the user’s manual for the unsteady transonic program used for this report. Reference 8 
describes the out-of-core direct solution subroutine. 



3.0 ABBREVIATIONS AND SYMBOLS 


f(x,y,t) 


l l 

ij,k 

I,J,K 


K 

le, LE 
M 

n,m 

q 

t 

TE, te 


X 0’y0 ,Z 0 


x,y,z 


Streamwise dimension of mesh region; also coordinate of downstream 
boundary 

Root semichord of wing or semichord of airfoil; also vertical dimension of 
mesh region 

Pressure coefficient, (p - p Q ) / (1/2 p 0 U Q ^) where p is the local pressure, 
p 0 the freestream static pressure, and p Q the freestream air density 

Instantaneous wing shape defined by z Q = 5f(x,y,t) 

Undisturbed wing or airfoil shape 

Unsteady contribution to wing or airfoil shape 

x,y,z subscripts and indices for points in the mesh 

AFT 

Transonic parameter, (1 - M^) / (M^e) 

Leading edge 

Freestream Mach number 

Mesh point indices 

co I e - ico (7 - 1 W, 

u xx 

Time in units of b / U Q ; also psuedo time defined by iterations in the complex 
differential equation for the unsteady potential 

Trailing edge 

Freestream velocity 

Physical coordinates, made dimensionless with the root semichord 

Scaled coordinates (x 0 ,py 0 ,pz 0 ) for the three-dimensional problem; the 

scaled coordinates for the two-dimensional problem are x and y, with x being 
the direction of fluid flow 



x ,y ,z 


x fie ,x te 


y 

P 



Ay?j 

4 % e 
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e 





£d?,f 




*1 

CJ 


Variables of integration 

Coordinates of leading and trailing edges 

TKy 

Vl -M 2 

Ratio of specific heats for air 

Jump in pressure coefficient 

Jump in <pj at plane of wing or vortex wake 

Jump in ipj , at wing trailing edge 

Thickness ratio or measure of camber and angle of attack 
(5 / M) 2 / 3 
coM / (1 -M 2 ) 

Critical value of Xj 

1/3 2/3 

Scale factor on y Q and z Q , p = 5 M ; also threshold pivoting factor (sec. 
6.1 and separation constant (sec. 7.0 and app. C) 

Coordinates for swept and tapered wing 

Complete, scaled perturbation velocity potential; also used for the unsteady 
potential in finite difference equations 

Steady scaled perturbation velocity potential 

Unsteady scaled perturbation velocity potential 

Angular reduced frequency (semichord times frequency in radians per second 
divided by the freestream velocity, cub / U) 


Matrix Notation 
[ ] Rectangular matrix 

{ } Column matrix 


[I] 


Unit matrix 



4.0 FORMULATION AND SOLUTION 


A detailed mathematical derivation of the method for the solution of the unsteady velocity 
potential for the flow about a harmonically oscillating wing is presented in reference 1 . The 
discussion here will be limited to a brief outline of the procedure for the two-dimensional 
flow. 


The complete nonlinear differential equation was simplified by assuming the flow to be a 
small perturbation from a uniform stream near the speed of sound. The resulting equation 
for unsteady flow is 

[k-(t- l)<P t -(7+ 1 )«p x ]«P xx + ^y y -(2<p xt + l ptt)/ e = 0 O) 

( 2 \ ( 2 \ 

where K = \1 “ M / / \ M e/, M is the freestream Mach number of velocity U Q in the 

x-direction, x and y are made dimensionless to the semichord b of the airfoil and the time 
t to the ratio b / U Q . With the airfoil shape as a function of time defined by the relation 

y 0 = §f(x,t) 

the linearized boundary condition becomes 

<Py = f x ( x ,t) + f t (x,t) ( 2 ) 

The quantity 6 is associated with properties of the airfoil (such as maximum thickness 
ratio, camber, or maximum angle of attack) and is assumed to be small. The coordinate 
y is scaled to the dimensionless physical coordinate y Q according to 

y = 5 1 / 3 M 2 / 3 y Q 

and e is given in terms of 5 by 

e = (5 / M) 2/3 

The pressure coefficient is found from the relation 

C p = - 2e(ip x + sp t ) 


The preceding differential equation is simplified by assuming harmonic motion and by 
assuming the velocity potential to be separable into a steady-state potential and a potential 
representing the unsteady effects. We write for a perturbation velocity potential 


s0 = s0cK x ’y) + ^l( x >y)e 


icot 


( 3 ) 


and for the body shape 

[ icutd 

f 0 (x) + fj(x)e J 

Since the steady-state terms must satisfy the boundary conditions and the differential 
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equation in the absence of oscillations, we obtain 


[K-(-, + lV Ox ]v> Oxx *v> Oyy = 0 


(4) 


with 


^0 y = f o( x ) . y = o - 1 < x< 1 (5) 

On the assumption that the oscillations are small and products of ipj may be neglected, 
equations (1) and (2) with the aid of equations (4) and (5) yield 

i|K -( 7 + 1 )<P 0 "|v?l l +¥>1 -( 2 ico/e)i p] +qip ]=0 ( 6 ) 

(L X J x ix yy x 

where 

2 , 

q ~ oj / e - ioo (7 - D^q xx 
subject to the wing boundary conditions 

=fi + icofi (x) , y = 0 -1<X<1 (7) 

y 

A computer program for solving the steady-state transonic flow about lifting airfoils based 
on equations (4) and (5) was developed by Krupp and Murman (refs. 9 and 10). The output 
of this program or a similar program can be used in computing the coefficients for the 
differential equation of the unsteady potential. The similarity of the unsteady differential 
equation to the steady-state equation suggests that the method of column relaxation used by 
Murman and Krupp for the nonlinear steady-state problem should be an effective way to 
solve equation ( 6 ) for the unsteady potential j . Note that equation ( 6 ) is of mixed type, 

being elliptic or hyperbolic whenever equation (4) is elliptic or hyperbolic. Central 
differencing was used at all points for the y derivative and all subsonic or elliptic points for 
the x derivatives. Backward (or upstream) differences were used for the x derivatives at all 
hyperbolic points. 

The boundary condition that the pressure be continuous across the wake from the trailing 
edge was found in terms of the jump in potential Aipj to be 

A n = A^ ]te e ' 10J ^ X ' Xte ) (8) 


where Ay> j is the jump in the potential at x = x te just downstream of the trailing edge and 

is determined to satisfy the Kutta condition that the jump in pressure vanish at the trailing 
edge. The quantity A<pi is also used in the difference formulation for the derivative i^i to 

yy 

satisfy continuity of normal flow across the trailing edge wake. 


For the set of difference equations to be determinate, the boundary conditions on the 
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outer edges of the mesh must be specified. In the original unsteady formulation, these 
boundary conditions were derived from asymptotic integral relations in a manner parallel to 
that used by Klunker (ref. 1 1) for steady flow. A later formulation (ref. 3) applies an out- 
going plane wave boundary condition to the outer edges of the mesh. This boundary 
condition is numerically simpler to apply and, on the basis of limited experience, appears to 
provide equally good correlation. 

The preferred numerical approach to solving the resulting large-order set of difference 
equations is a relaxation procedure, which permits the calculation to be made as a sequence 
of relatively small problems. However, as discussed in preceding NASA reports by the 
authors (refs. 2 and 3), a significant problem of convergence with the relaxation pro- 
cedure was encountered that severely limits the range of Mach number and reduced 
frequency for which solutions may be obtained. The authors currently believe the only 
practical technique for circumventing these instabilities is a full direct solution where the 
difference equations are solved “all at once” rather than by line relaxation. 
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5.0 REPRESENTATION OF THE PRESSURE PULSE 


Consider an airfoil submerged in a flow that is subsonic but of relatively high Mach number. 
With appropriate conditions, the usual subsonic pressure distribution differs from that at 
lower Mach numbers by the presence of shocks attached to the airfoil. For the simple case 
of a symmetrical airfoil, and with the airfoil at zero angle of attack, these shocks would be 
of equal strength and symmetrically placed top and bottom. As the airfoil oscillates in 
pitch, the shocks move in opposite directions, resulting in a pressure pulse that will oscillate 
in amplitude. In this particular case, although the shocks move fore and aft, the resulting 
pressure pulse is fixed at a point on the chord. It is also possible, of course, for the location 
of the pulse to move if the airfoil section is not symmetrical or if the airfoil is at a nominal 
angle of attack. Also, if the oscillation is very slow, the pulse is in phase with the motion. 
As the frequency is increased, oscillation of the pulse may lag the motion of the wing. This 
pressure pulse, called a “pressure doublet” by Ashley, is described in some detail by 
Tijdeman (ref. 4) and analyzed in a general fashion for its potential effects on flutter by 
Ashley (ref. 5). The important characteristics of the pressure pulse with respect to flutter 
instabilities appear to be its amplitude, its phasing with respect to the section motion, and 
its location relative to the section elastic axis. In particular, for airfoils with blunt leading 
edges such as supercritical sections, the pulse tends to be forward on the section, which 
increases the likelihood of instabilities. 

An illustration of how the moving shock pulse may be represented in the frequency domain 
solution of the type discussed here is obtained by considering the Fourier expansion of the 

pressure pulse produced by two unit step functions oscillating back and forth 90° out of 
phase on opposite sides of a plate. The resulting lift pulse oscillates in width and the center 
position oscillates along the plate. In addition, the amplitude of the pulse changes sign in 
the manner of a square wave. These properties are easily seen in the graphs of figure 4. The 
Fourier coefficients of the sine and cosine of the first harmonic are shown as a function of 
position along the plate in figure 5. This is seen to be very similar to the sharp peaks 
observed in figure 2 which appear at the mean shock location in figure 3. Thus it appears 
that the present frequency domain method may well produce the pulse amplitude and 
movement of the fundamental resulting from the moving shock system without additional 
treatment such as computing the actual shock motion by a nonlinear theory and performing 
the Fourier analysis. The derivation of the equations for this analysis is included in 
appendix D. 

Tijdeman (ref. 4) has performed a similar analysis for a single shock wave on one side of the 
airfoil. He showed that the major part of the unsteady pressure pulse is due to the first 
harmonic. 





6.0 THE DIRECT SOLUTION 


The full direct solution of the finite difference formulation of equation (6) was explored in 
reference 3 as a means for circumventing the solution instability problems associated with 
the relaxation procedure. These initial studies were successful in that solutions could be 

obtained for values of X j = ocM / (l - M^) well above the critical value. However, the accu- 
racy of the solutions (as measured by comparing solutions for a flat plate with subsonic kernel 
function solutions) rapidly decreased with increasing X j . This accuracy problem was at first 

assumed to be due to an insufficient number of mesh points resulting from the use of an 
in-core solution routine. An out-of-core solution has been developed under the current 
contract, and the number of mesh points to be handled can now be increased to the number 
used in the relaxation procedure. 

The new program has proved successful although all the improvement in correlation has not 
been due solely to increasing the number of mesh points. Other modifications that have 
contributed significantly to obtaining better results include removing a coordinate trans- 
formation, which provided for applying the boundary conditions at infinity to the outer 
boundaries of the finite difference region, and redistributing the mesh points vertically to 
more uniformly span the solution region. 

The out-of-core solution procedure, which is developed in detail by Yip (ref. 8), is discussed 
in the following section. 


6.1 THE OUT-OF-CORE SUBROUTINE 

The out-of-core solution routine is set up to solve the matrix expression 

[A]{x} = {R} (8) 

The coefficient matrix, [A] , is complex, of the order 2000 to 4000, sparse and banded. 

The right-hand-side matrix, {R}, is a column matrix that introduces the internal boundary 
condition (i.e., the boundary conditions on the wing) and for the general flutter problem 
may have as many as 20 different values for each [A] . The forms of the matrices 
of equation (8) that result from using outgoing wave boundary conditions are shown in 
appendix A. 

The coefficient matrix, [A] , of equation (8) can be written as the sum of two matrices 

[A]=[B]+[D] (9) 

where [B] is a banded matrix and [D] is such that its upper triangular part is all zero, and 
its lower triangular part has four columns of nonzeroes. The matrix [D] can be written as 

[D] = [U] [V] T (10) 
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where [U] and [V] are nx4 matrices. If the solution is for 


[A] { x } = {R} (11) 

then, by a variant of the Sherman-Morrison formula, we have 

{x}=[A]' 1 {R} = [B]" 1 {R} - [B ] _1 [U] [l 4 +V T B _ 1 u] 1 [V] T [B]' 1 {R) (12) 

T T-l I 

upon noting that II4 + V B UJ is a 4x4 matrix. 

A two-step method is used where 


1. [B] is decomposed into a product of a lower and an upper triangular matrix 
by an algorithm for block banded (or profile) matrices. Then 

[B] { R } and [B] ^ {U } are solved by forward and backward substitutions. 

2. Equation (12) is used to compute [A] {R}. 


The matrix is partitioned into blocks (as described in app. A) with the order of each block 
equal to the number of vertical mesh points. Only nonzero blocks are stored as a record in 
a random access file. Since any algorithm for LU factorization of banded matrices ultimate- 
ly fills in zeroes between bands, all the nonzero blocks of the coefficient matrix are stored 
as full matrices. A “threshold” pivoting strategy is used. That is, a nonnegative number, p, 
is chosen so that 

0<p< 1 


and if 


| a kk[ • | a jk| J >k 

then there is no pivoting in the k tk pivotal step. If p is set equal to zero, there is no pivoting; 
if p = 1 , there is the usual row pivoting. Our experience shows that using p = 0.01 has been 
very satisfactory. 


The numerical stability of the solution procedure is monitored by comparing the actual 
solution with a computed solution. A new right-hand side { R} can be generated where the 
element Rj is the row sum of the i tk row of [A] . Therefore, if {x }is the actual solution of 

[A] {x} = { R }, then the deviation of { X } from { 1 } will give some indication of the quality 

-10 

of the solution. In all solutions run so far, this error has been of the order 10 . The 

program has the capability of treating several wing motions for the same frequency and Mach 
number with considerable economy in computing resources. 


6.2 PROCEDURAL MODIFICATIONS 

As will be noticed by examining section 6.3, the direct solution program provides relatively 
good solutions for the two-dimensional airfoil section of vanishing thickness at relatively 


large values of Xj, i.e., for values well above the critical value of Aj for relaxation procedures. 

This represents a significant improvement over the direct solution results discussed in 
reference 3. Although development of the out-of-core solution subroutine is the largest 
single modification made during the current work, several other changes have contributed to 
the better results. 


First, the coordinate transformation that applied the physical boundary conditions at infinity 
to the outer boundaries of the finite difference solution region failed to yield good results 
for higher A j values and was removed. Although this transformation — which works well for 

steady state and small values of Aj but does not work for larger values of Aj — has not been 

studied, we assume that the problem lies with trying to analyze a set of waves extending to 
infinity with a limited number of mesh points. That is, it appears to be important to have 
enough mesh points to define properly all the waves within the finite difference region. 

This is demonstrated by the better results obtained by spreading the mesh points out rather 
uniformly in the crossflow direction. Points are still stacked close to the wing surface where 
the gradients in <pj are known to be large. However, increasing adjacent intervals going away 
from the wing by a factor of 1 .4 did not provide a good point distribution. 


Outgoing wave boundary conditions have been used and appear to be adequate. However, 
the equations given in reference 3 for the upper and lower boundaries are in error and 
should be as follows: 


icoMYK 

^1 + 2 ~ *1 =0 

y 1 -M 


UPPER BOUNDARY 


icoMVK 

V ?1 - 2 " *1 -0 

y 1 -M 


LOWER BOUNDARY 


An alternative downstream boundary condition has been developed and applied. This 
involves replacing the outgoing wave boundary condition with the value ipj on the down- 
stream boundary as determined from the flow induced by the doublet sheet shed from the 
trailing edge of the airfoil. The calculated pressure distributions at M =0.9 and co = 0.6, for 
example, are essentially the same for both downstream boundary conditions. 


There have been no apparent problems with numerical instabilities resulting from the 
presence of eigenvalues of the coefficient matrix. The pressure distributions appear to be 
smoothly varying with small changes in reduced frequency. 


6.3 ANALYTICAL RESULTS 


This section presents and discusses some results from the direct solution program described 
in the preceding section. These examples include both subsonic and supersonic freestream 
conditions and are for airfoil sections only. That is, all the examples are two-dimensional. 
The following section examines the purely subsonic problem in some detail as it is impera- 
tive that good correlation be obtained here if the full mixed flow condition is to be 
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successfully handled. The succeeding section presents some supersonic freestream results 
for a circular-arc airfoil. 

6.3.1 SUBSONIC FREESTREAM EXAMPLES 

The flow around an airfoil of vanishing thickness is purely subsonic and corresponds to the 
problem solved by conventional subsonic unsteady flow programs such as that developed 
for NASA by Rowe et al. in references 12 and 13. This latter program is considered to 
provide very accurate pressure distributions, and the finite difference results should 
approach those of the current NASA study as the number of mesh points increases to 
infinity. Thus, the criterion for evaluating pressure distributions is the correlation of the 
finite difference calculations with the results from the NASA program of references 12 and 
13. 

The first point to be made is that good pressure distributions have been obtained from the 
direct solutions for values of X j well above the critical value for which relaxation solutions 
fail to converge. Results for a section oscillating in pitch about the leading edge at a free- 
stream Mach number of 0.9 and four reduced frequencies are presented in figures 6 through 
9. At this Mach number and for our standard sized finite difference solution region, the 
critical reduced frequency would be about 0.1 . As can be seen, very good correlation is 
obtained up to a reduced frequency of 0.9. 1 


Results for a second example for the same motion at Mach 0.4 are shown in figures 10 and 
1 1 . Here, correlation with the more exact theory is good at both co = 0.9 and cu = 3.0 The 
critical value of Xj is about 1.0. 

Finally, results are presented for a section having an oscillating quarter-chord control surface 
at M = 0.5 in figures 12, 13, and 14. Here, correlation is seen to be good for reduced fre- 
qencies of 0.5, 1 .0, and 1.6. 

Experience to date indicates a number of characteristics are important in setting up a 
satisfactory mesh point distribution. Nonuniform mesh point distributions are used in order 
to concentrate points where the gradients in i are largest and still keep the number of points 
to a minimum for calculation purposes. In the x-direction (direction of flow), a high 
density of points are located in the vicinity of the airfoil leading edge, the control surface 
hinge-line, and the trailing edge where the Kutta condition is to be satisfied. Having met 
these conditions, the pressure seems relatively insensitive to the particular distribution of 
points. In the crossflow direction, grid points are located very close to the wing surface to 
accurately define y?j directly adjacent to the wing surface. It is also important to keep enough 
points distributed in the outer parts of the flow to properly define the wave shapes. 

An example of a distribution with good correlation at relatively high values of Xj is shown 
in figure 15. The complete grid is not shown, but it is symmetric about a line halfway 
between the j = 30 line and the j = 3 1 line. Here, the finite difference points are the inter- 
sections of the lines. The yjj distributions for this grid with the airfoil at a Mach number of 
0.9 and four different reduced frequencies are shown in figures 16 through 19. A descrip- 
tion of the ipj plots is given in appendix B. The potential, which in this case is antisymmetric 
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Figure 8. -Jump in Pressure Coefficient Across a Flat Plate Oscillating in Pitch, 
M = 0.9, co = 0.6 





Figure 9.— Jump in Pressure Coefficient Across a Fiat Plate Oscillating in Pitch; 
M = 0.9, u) = 0.9 
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Linear theory 
(NASA CR-2543) 



Figure 12 — Jump in Pressure Coefficient Across a Fiat Plate With Oscillating 
Quarter-Chord Control Surface; M = 0.5, gj = 0.5 



Linear theory 
(NASA CR-2543) 



Figure 13 — Jump in Pressure Coefficient Across a Flat Plate With Oscillating 
Quarter-Chord Control Surface; M = 0.5, co = 1.0 



_ — _ Linear theory 

(NASA CR-2543) 



Figure 14.— Jump in Pressure Coefficient Across a Flat Plate With Oscillating 
Quarter-Chord Control Surface; M = 0.5, oj = 1.6 
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J=40 



1=1 1=10 1=37 1=56 


Leading edge between X(l = 10) and X( 1 = 1 1 ) 
Trailing edge at X(l=37) 

Airfoil between Y(J=30) and Y(J=31) 


Figure 15. — Mesh Pattern for 56 x 60 Grid 


29 






Figure 18.- 
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Figure 19.— Velocity Potential for a Fiat Plate Oscillating in Pitch; M = 0.9, co = 0.9 
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about y = 0, is continuous in the crossflow direction ahead of the leading edge. The discon- 
tinuity in y>j across the airfoil from which the lift is calculated shows clearly, as does the 
continuing discontinuity across the wake. The boundary conditions across the wake ensure 
that the pressure difference across the wake is zero. 

It will be noted that the gradients in the flow direction, even adjacent to the plane of the 
wing (there are no mesh points along y = 0), are, in all cases, relatively mild except at the 
airfoil leading edge. The most noticeable characteristic of the distribution is the develop- 
ment of bow-type waves, with the wavelength decreasing as the reduced frequency is in- 
creased. The waves appear to be well defined by the grid of figure 15 for co = 0.3 and co = 0.45. 
0.45. However, waves at co = 0.6 are not as smooth as those obtained for smaller frequencies. 

Several different grids were used for co = 0.45 and co = 0.6 in an attempt to improve the 
pressure distribution results. The pressure results were relatively insensitive to movement 
of the upstream and downstream boundaries, as well as to the addition of points close to 
the upper and lower boundaries. 

A grid resulting in less accuracy is shown in figure 20 for which rows of mesh points are 
stacked very closely against the wing and spread out toward the upper and lower boundaries. 
Although good results are obtained at small values (e.g., co’= 0.18 at M = 0.9), the quality of 
results decreases rapidly with increasing reduced frequency. Also, at the higher values of 
reduced frequency where the pressures are bad, the results are sensitive to the location of 
the downstream boundary. The velocity potential for the grid of figure 20 and co = 0.9 and 
co = 0.3 is shown in figure 21 . Although the velocity potential is smooth in the region 
adjacent to the airfoil, the wave pattern in the outer field is not well defined and is distinct- 
ly different from that of figure 16. The downstream flow pattern appears different because 
the downstream boundaries are at different locations; 4.0 semichords aft of the trailing edge 
in figure 16 and 1.75 semichords in figure 21. A great variety of point distributions in the 
x-direction were tried with this type of distribution in the crossflow direction. The resulting 
pressures proved to be sensitive to the location of the downstream boundary. Also, good 
correlation was never achieved for this combination of Mach number and reduced frequency 
using this crossflow pattern. 

The velocity potential for another pitch case is shown in figure 22 with the grid of figure 15. 
Here the Mach number is 0.4 and the reduced frequency is 3.0. The X j for this case matches 
that for M = 0.9 and co = 0.3. The upper and lower boundaries are at the same scaled distance 
from the airfoil surface (6.25 semichords above and below the wing), but this length transforms 
into a larger physical length at M = 0.4 than at M = 0.9: hence, the additional number of waves 
in the crossflow direction. Despite the number of waves, the pressure distribution correlates 
well with more exact theory (see fig. 1 1). 

The velocity potential for a section of vanishing thickness with an oscillating quarter-chord 
control surface is shown in figure 23. The freestream Mach number is 0.5 and the reduced 
frequency is 1 .6. The downstream boundary is 9 semichords aft of the trailing edge. Again, 
despite the number of waves in the crossflow direction, correlation with more exact theory 
is good, (see fig. 14). 
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1=1 1=19 1=48 1=64 


Leading edge between X(l = 19) and X{ 1=20) 
Trailing edge at X ( 1=48) 

Airfoil between Y(J=25) and Y(J=26) 


Vertical spacing: Y(J=27)-Y(J=26)=Y(J=26)-Y(J=25)=Y(J=25)-Y(J=24) 
Y(J)-Y(J-1)=1.4«[Y(J+1)-Y(J)] 2<J<24 

Y(J+1 ) Y(J)=1 ,4*[Y(J)-Y(J-1 )] 27<J<49 


Figure 20 — Mesh Pattern for 64 x 50 Grid 
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6.3.2 SUPERSONIC FREESTREAM EXAMPLES 


The pilot two-dimensional program with the full direct solution procedure has been adapted 
to the supersonic freestream condition. For demonstration purposes, sample calculations 
have been made at two Mach numbers and for a range of reduced frequencies. The results 
appear to be reasonable, although we have nothing with which to correlate them. No 
numerical problems were encountered. 

As a trial example, pressure distributions were calculated for a 6% thick circular-arc airfoil 
oscillating in pitch with freestream Mach numbers of 1 .05 and 1.15 and for a range of 
reduced frequencies up to 0.72. The steady-state pressure and velocity potential distribu- 
tions were calculated using a Boeing routine that was derived from a program developed by 
Murman at NASA Ames (ref. 14). The procedure is fully conservative. Plots of pressure 
coefficients for the M = 1.15 calculations are given in figures 24 through 27. 

The steady-state pressure coefficient is given in figure 24. The corresponding data from 
Traci et al. (ref. 15) is also shown, with the difference in location of the bow shock between 
the two procedures being attributed to differences in point spacing ahead of the leading 
edge, although this effect has not been investigated. 

The unsteady pressure distributions for co = 0.06, 0.18, 0.36, and 0.72 are presented in figure 
25. The results include the singularity in the pressure distribution at the leading edge, 
reflecting the effect of the subsonic bubble attached to the leading edge. The pressure aft 
is relatively constant as would be expected in purely supersonic flow. 

The steady-state velocity potential is presented in figure 26. The size of the finite difference 
region used for the unsteady calculation is much larger than that used for the steady calcula- 
tion to ensure that the boundary conditions used on the outer boundary did not affect the 
pressures on the wing. The region used for the steady calculation is small enough to 
truncate the full crossflow extent of the i distribution. This truncation occurs only in the 
regions that do not influence the pressures on the wing itself. 

The unsteady velocity potential distribution for M = 1.15 and co = 0.06 is shown in figure 27. 
This distribution is relatively smooth considering the small number of grid points in the 
portion of the flow that actually affects the airfoil section. There are a couple of sharp 
breaks in the distribution at the outer edges of the region of influence. One irregularity 
shows clearly in the plot of the the imaginary part. It is assumed these breaks are due to the 
truncation of the ip Q distribution discussed above and the relatively small number of mesh 
points used. 

The results for the M = 1 .05 calculation are presented in figures 28 and 29. A plot of the 
pressure coefficient for steady flow is shown in figure 28. The subsonic bubble about the 
leading edge extends further in the upstream direction than that for the M = 1.15 case, but 
only slightly further aft of the leading edge. The unsteady pressures are shown in figure 29 
for co = 0.06 and 0.12, and have the same characteristics as those at M - 1.15. 
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Without shock point differencing 
4- — With shock point differencing 
O Boeing solution 


Data from ref. 1 5 





GO 


Thickness ratio is 6% 


Figure 24— Jump in Pressure Coefficient Across a Circular-Arc Airfoil; M = 1.15 
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co = 0.18 
03 = 0.36 
• • • • co = 0.72 

a. REAL PART 

Figure 25 — Jump in Pressure Coefficient Across a Circular-Arc Airfoil Oscillating in 
Pitch; M= 1.15 




05 = 0.06 

•— • 05 = 0.18 

• 05 = 0.36 
• • • • 05 = 0.72 


b. IMAGINARY PART 


Figure 25 — (Concluded) 
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Figure 26 — Steady-State Velocity Potential for a Circular-Arc Airfoil; M = 1.15 
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Figure 28.— Pressure Coefficient Distribution for a Circular-Arc Airfoil; M = 1.05 
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Figure 29 — Jump in Pressure Coefficient Distribution for a Circular-Arc Airfoil Oscillating in 
Pitch; M= 1.05 
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6.4 FUTURE APPLICATIONS 


The frequency limitation on convergence occurs for the relaxation solution of three- 
dimensional transonic unsteady harmonic flow as well as two-dimensional flow, although it 
is not quite as severe due to the extra dimension. Nevertheless, to cover the required range 
of frequency, a direct solution seems to be the only practical approach. 

Approximate solutions of the three-dimensional problem may be obtained by using a strip 
theory approach where the two-dimensional direct solution is applied to each spanwise 
cross-section, with the steady-state potential derived from the actual three-dimensional 
steady-state transonic small-perturbation method of Bailey and Ballhaus (ref. 16). 
Reference 6 discusses this procedure in greater detail regarding its application to aeroelastic 
analyses. 

The matrix derived from the system of difference equations for the full three-dimensional 
problem is extremely large. For example, if the numbers of mesh points in the x, y, z 
direction are, respectively, 50, 40, 60, then the order of the coefficient matrix is 120 000. 

If the present out-of-core linear equation solver is used, the immediate problem will be 

10 

insufficient disk storage. In the present example, it will take 1 .06 x 10 sectors to store 
the factored matrix. Another problem would be the amount of central processor unit 
(CPU) time required for the computations. For this example, it would take approximately 
5 

144 x 10 CPU seconds on the Cyber 175 (assuming one multiplication takes 3.4 nano- 
seconds). 

However, there are several possible practical approaches to solving the matrix equation for 
the three-dimensional problem. It should be noted that the generalized conjugate gradient 
method provides a means for solving large matrix equations with practically no disk storage 
and within a reasonable amount of CPU time. The generalized conjugate gradient method, 
like iterative procedures, uses the original matrix in such a way as to eliminate the problem 
of “fill-ins” of zero elements, as in the case of Gaussian elimination. This is a very good 
feature for the application under discussion because the original matrix can be recomputed 
from time to time quite economically. In the example cited here, it is estimated that it 
would take 576 CPU seconds to arrive at the solution (excluding the CPU time taken to 
recompute the matrix). However, the biggest concern in using the generalized conjugate 

T 

method is that the condition number of the required matrix [A] [A] is the square of the 
condition number of [A] . This means that if [A] is only fairly well conditioned, the 
solution obtained might not be satisfactory. Yet, recently, numerical analysts have been 
studying preconditioning techniques (e.g., ref. 17) to improve the condition number of [A] 
before applying the generalized conjugate gradient method. 

The vector processors such as the CDC Star 100 and the more advanced CRAY machine 
may also provide a practical means for obtaining direct solutions to the unsteady transonic 
problem. The difficulty with the application of vector processors to the matrix of the 
difference equations is the sparseness of the resulting nonzero vectors of the coefficients. 
These machines are especially suited for large systems of equations for which the matrices 
are full, and they may be the only feasible approach to such large-scale problems. Rodrique, 
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Madsen, and Karush (ref. 18) propose an odd-even reduction method for banded linear 
equations. The three-dimensional problem, unlike the two-dimensional one, has a strictly 
banded matrix. If the x, y, z grid is 40 x 40 x 60, then the band width is 50 x 60 x 3 = 9000 
and the method of Rodrigue et al. is applicable. By the application of permutation matrices 
P, the odd variables are decoupled from the even variables. For example, consider an even- 
order matrix equation 

[A] { x } = {b } 

T 

If [P] is a permutation matrix and [P] [A] [P] has the form 

[P] [A] [P] T = 

then the approach is to construct a matrix [Q] such that 

[Q] [P] [A] [P] T = 

where [D] is an easily invertible matrix and [A] has the same band structure as [A] . The 
solution of [A] {x} = {b} is equivalent to 

[Q] [P] [A] [P] T [P] {X} = [Q] 




Letting { y } = [P] { x } 


y l 

y 2 


and {b}= [Q] [P] <b } 



the solution of the system 

[Q] [P] [A] [P] T { y > = (b) 

is seen to be 

jy 2 | = [A] 1 jb 2 | 

)y 1 j-[D] _, (jb 1 [-[Uljy 2 j) 
and the final solution of[A]{x} ={b} is given by 

{x} = [P] T {y} 

The system [A] { Y} = jb 2 j still needs to be solved, but it is of the same form as the original 

equations and hence the odd-even reduction may be applied to this set to reduce still further 
the order of the system. Reference 18 presents some sufficient conditions which guarantee 
that odd-even reduction can be performed, but conditions which guarantee that the algo- 
rithm can be applied in a cyclical manner still need to be established. The overall stability 
of the method, moreover, needs to be analyzed theoretically. 
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7.0 AN EXACT TWO-DIMENSIONAL ANALYSIS 


In reference 3, the one-dimensional version of the flat plate problem for unsteady transonic 
flow was solved by a finite difference procedure and the results compared with the analyti- 
cal solution. One of the main results of this study was that the outer mesh boundary condi- 
tions could significantly affect the accuracy of finite difference solutions. This result may 
be described as follows: The unknown constants in the general solution for the exact 
partial differential equations are determined by using the boundary conditions to write a set 
of simultaneous equations. This set of equations forms a characteristic value problem with 
Aj as an eigenvalue. The characteristic values of Aj may be real and/or complex. In the 
finite difference procedure, solutions are calculated only for real values of A j . A plot of error 
(the maximum difference between the analytical and finite difference solutions) versus oo 
is shown in figure 30, which is reproduced from reference 3. It is seen that if the boundary 
conditions are such that the eigenvalues (the characteristic values of A j) are complex, a 
smooth error curve is generated. If the boundary condition results in eigenvalues that are 
real, the error curve tends to infinity at values of co corresponding to the eigenvalues. A 
parallel study for this report shows that this same phenomenon also occurs for the two- 
dimensional problem. 


The general form of the equation for a flat plate oscillating in two-dimensional flow is 


Kipi + 1^1 -2i(oj/e)v?i + (cu / e)ipi = 0 
*xx *yy x 1 


(13) 


Three exact solutions of equation (13) are derived in appendix C for the purpose of study- 
ing the effects of mesh boundary conditions on the accuracy of finite difference solutions. 
These are : 

<Pj = jexp[iA]M(x - xj)] • sin(7ry / b) • sinh[q(a - x)] [/ sinh(q(a - Xj)J, Aj <n / (bVK) 
= jexp [iAjM(x - xj)] • sin(7ry / b) • sin[q(a - x)] (/ sin [q(a - Xj)J , Aj > n / (bVK) 

q cosh[q(a - x)] + iA] sinh[q(a - x)] ) 


= exp [iAjM(x - Xj) ] • sin(7ry / b) •< 
= exp[iAjM(x - xj)] • sin(7ry / b) 


'q cosh [q(a - Xj)] + iAj sinh[q(a - x j)] ) 
q cos[q(a - x)] + iAj sin[q(a - x)] ) 

q cos [q(a - x j)] + iA j sin [q (a - x j ’ 


where q = ~\]\^\^ -i? I (b^K)| 

i pj = exp[iAjM(x - x j)] • £exp (iA j zy Vk) + ^ exp^- iAjzy VK^j 

((Ml - Aj) exp[ipj(x - a)] + (qj + Aj)exp [- ipj(x - a)] ) 


Aj < 7T / (bVK) 

Aj >7T/(bVK) 

(14) 

(15) 


((q t - Aj)exp[iqj (xj - a)] +(qj + Aj)exp[- ipj (xj - a)] 


(16) 
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Figure 30.— Comparison of Error Curves for Boundary Conditions Yielding All-Real and 
Complex Eigenvalues 
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( 17 ) 


and /ij = iA 



and z is a root of 



The coordinate x = xj is the upstream boundary, x = a the downstream boundary, while 
y = 0 and -b are the upper and lower mesh boundaries. 

Three types of boundary conditions were considered. In all three cases, the unsteady poten- 
tial was prescribed on the upstream boundary. For the first solution (eq. (14)), the pertur- 
bation potential was set to zero on the other three boundaries. For the second solution 
(eq. (15)), the potential was set to zero on the upper and lower boundaries, and outgoing 
wave boundary conditions were applied on the downstream boundary. For the third solu- 
tion (eq. (16)), outgoing wave boundary conditions were prescribed on all three boundaries. 
The first leads to real eigenvalues, the second leads to imaginary eigenvalues, while the third 
set of boundary conditions yields complex eigenvalues. 

The direct finite difference solutions are compared with the exact analytic solutions in 
figure 31. As expected, for the first case the error peaks at those reduced frequencies that 
correspond to the eigenvalues of the zero potential boundary conditions. Surprisingly, con- 
siderable error also occurs for the second solution about these same resonant frequencies, 
although the eigenvalues for this solution are imaginary. 

Outgoing waves on all three boundaries eliminate this difficulty, and the error in the third 
solution is smooth and considerably less than that for the solution with the zero potential 
boundary conditions. The improved accuracy for the outgoing wave boundary conditions 
on all three boundaries is even better than the graph indicates. The curves in figure 3 1 are 
absolute magnitude of error. Since the maximum value on the upstream boundary for the 
first two solutions is unity while the maximum value for the third solution is about three, 
all three boundary conditions yield about the same relative error in the range of reduced 
frequency: 0.1 to 0.2. 
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Figure 31.— Influence of Mesh Boundary Conditions on the Accuracy of the Finite 
Difference Solution of the Linearized Differential Equation for 
Unsteady Subsonic Flow 
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8.0 ADDITIONAL STUDIES 


Additional studies performed as part of the investigations for this report are presented in a 

separate report (ref. 6). These studies are either theoretical in nature and were not imple- 
mented, or are concerned with concepts that were tried and did not significantly affect 

either the solution or the solution procedure. 

1. An analysis is made of the finite difference operator applied to supersonic mesh 
points. The usual five-point upwind difference operator is found to be over-stable. 
This means that, although the numerical solution approaches the true solution as the 
mesh (step size) is refined, for a given mesh the correct numerical solution is distorted 
by being attenuated as the disturbance propagates through the mesh. A new nine- 
point operator is derived which is stable and thus does not attenuate the disturbances. 

2. The effect on relaxation solution convergence of adding a viscous term is investigated. 
It is found that the addition of viscosity has very little effect on convergence, even 
when large amounts of viscosity are used. 

3. An alternate downstream boundary condition is applied in which the velocity poten- 
tial on the downstream boundary is set equal to the potential due to the wake. The 
resulting pressures are not significantly different from those obtained by assuming a 
plane outgoing wave. 

4. The three-dimensional differential equation for unsteady transonic flow, as used in 
this report, is rederived using an oblique coordinate system appropriate for use with 
swept and tapered wings. 

5. In reference 2 it was shown that for the two-dimensional problem, row relaxation 
converged more rapidly than column relaxation, but that additional terms were 
required for mesh points at which the flow equation was hyperbolic. These addition- 
al terms are derived for the three-dimensional problem. 

6. Equations are derived and presented for implementing a procedure that matches the 
finite difference solution over an inner mesh region to an outer linearized solution 
that has the appropriate outgoing wave properties. This procedure should reduce the 
size of the finite difference region as well as provide a more accurate mesh boundary 
condition. 

7. A three-dimensional aeroelastic solution procedure is described that uses two- 
dimensional unsteady transonic air forces at each cross section. The two-dimensional 
air forces are calculated using the three-dimensional steady-state velocity potential 
distribution. 
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9.0 CONCLUSIONS 


The direct finite difference method of this report appears to yield sufficient accuracy for 
two-dimensional subsonic flow at Mach numbers and reduced frequencies of interest to 
flutter analysts. The most pressing problem at the start of the investigation had been lack of 
accuracy at higher values of the parameter Xj. As shown by examples in section 6.0, this 
inaccuracy seems to have been overcome by removing a coordinate transformation that 
placed the outer boundaries of the finite difference region at infinity, by using a more uni- 
formly spaced grid in the crossflow direction, and by increasing the number of finite diff- 
erence points to provide a finer grid in the field away from the wing. The correlation of the 
finite difference method with the more exact linear theory is good and no numerical diffi- 
culties have been encountered. As a result of this investigation, we have also concluded that 
the simple outgoing plane wave boundary condition yields sufficient accuracy. 

Two concerns of immediate interest remain. First, there is a question of correlating the 
analytic results of the harmonic procedure with valid experimental data. Preliminary corre- 
lation studies with the experimental studies of Tijdeman (ref. 4) have been encouraging but 
not conclusive. The forthcoming availability of the NASA-Ames data for two-dimensional 
sections provides an opportunity to further these studies. Second, there is the problem of 
extending the direct solution to three-dimensional flow. The two possible procedures dis- 
cussed in section 6 warrant further investigation. 

There is controversy over whether the harmonic procedure as formulated for this report 
includes the characteristics of transonic flow that are of most importance to flutter instabil- 
ity. As of now, the main concern is whether the effects of moving shocks are properly in- 
cluded in the calculation. The most obvious effect — the oscillating pressure pulse that 
results from the moving shock system — is contained in the numerical results. We believe 
that this pulse is the first harmonic of the Fourier expansion of the solution for the full 
nonlinear system and probably incorporates the most significant effects of the moving shock 
system with respect of the flutter instability. This is an important point and certainly 
worthy of further investigation. 
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APPENDIX A 


FORM OF DIRECT SOLUTION MATRICES 


The matrix of coefficients, [A] , of equation (8), is a large, relatively sparse, banded matrix 
with complex elements. The two-dimensional solution region with the index conventions 
used in reference 1 is shown in figure A-l . Each mesh point interior to the boundaries is 
numbered in sequence starting with the first point inside the lower left-hand (upstream) 
corner. The points on the boundary are eliminated with the outer boundary conditions. 
The mesh points are numbered by column, bottom to top, and upstream to downstream. 
The resulting coefficient matrix is of the order* (IMAX - 2) • (JMAX - 2) and is conven- 
iently partitioned into blocks, with each block associated with a column of mesh points 
from the solution plane. There are thus (IMAX - 2) • (IMAX - 2) blocks each of order 
(JMAX - 2). The block matrix is shown in figure A-2. Block matrix [A] is a tridiagonal 

matrix with elements that are coefficients of <^i , , and <p-\ . Arrangement of the 

i,j-l ij MJ+1 

elements in [A] is shown in figure A-3. Block matrix [B] is a diagonal matrix whose non- 
zero elements are coefficients of <p i . Block matrix [C] is a diagonal matrix whose 

i-I,j 

elements are coefficients of ipi . If a mesh point is supersonic, as determined from the 

i+l, J 

steady velocity potential distribution, the corresponding diagonal element in [C] is zero. 

The block matrix [D] is a diagonal whose elements are coefficients of ipi . If a mesh 

i-2,j 

point is subsonic, the corresponding diagonal element in [D] is zero. Block matrix [E] , 
illustrated in figure A-4, is a diagonal matrix with elements from the wake calculation in- 
cluded. Block matrix [F] has four nonzero elements as shown in figure A-5. Block matrix 
[G] is shown in figure A-6. 


* The dimensions given in this paragraph are for full-space solutions. For half-space 
solutions, (JMAX - 2) is replaced by (JM -1). 
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Figure A-3.— Form of Block Matrix A 
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Figure A-4.—Form of Block Matrix E 
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Figure A-5— Form of Block Matrix F 
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Figure A-6.— Form of Block Matrix G 
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APPENDIX B 


PLOTS OF VELOCITY POTENTIAL DISTRIBUTIONS 


A capability for plotting the complex velocity potential distributions calculated with the 
pilot two-dimensional program of this report has been developed. This capability has been 
implemented on the PDP-1 1/70 computer system. These plots have permitted visual study 
of the potential distributions and provide a basis for estimating the number of points re- 
quired and the distribution of the points within the solution space for a given case. The real 
and imaginary parts of the velocity potential are plotted as separate distributions. Each 
distribution is plotted as a three-dimensional surface over the solution region. The figures 
appear on a scope and may be manipulated by the observer (i.e., rotated about any of the 
three axes and regions of particular interest may be enlarged) in order to examine the 
characteristics of the distribution. In practice, the distributions for the problems of this 
study are so complicated that no single view provides a clear understanding of the total 
field, and the ability to change the view continuously has proved invaluable. Hard copies 
may be obtained ; examples are included in section 6.0. A schematic of these plots is shown 
in figure B-l. 
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Figure B-1.~ Schematic of Velocity Potential Figures 


APPENDIX C 


EIGENFUNCTION ANALYSIS FOR A FLAT PLATE 
IN TWO-DIMENSIONAL FLOW 

The partial differential equation for the harmonically oscillating flat plate (see eq. (6)) is 

- 2i(c o / e) i£] + (oj /e^y?]=0 

Substituting K ■= (l - M )/M~e yields 


(C-l) 


(C-2) 


2 2 2 

2iM u w M 

<Pi — v?i + — V i -0 

*xx ‘yy 2 x 2 

1 -M 1 -M 

where y = Yk” y. We can eliminate the first derivative term by substituting 

f / 2 2\ 1 iXiMx 

= expl_i\6jM /|3/x Ji// = e \jj 

2 2 2 

where X] = coM / p and p = 1 - M . Substituting equation (C-2) into equation (C-l) 
reduces the differential equation to 

^xx + ^yy + ^l ^ = 0 (C’ 3 ) 

The boundary condition representing an outgoing plane wave in the x-direction is given by 

= 0 (C-4) 


i^i + ico 
x 


or 


i P] + iXj(l - M)ipj = 0 


Since = ^iXjMi// + i// x Je^^ X , the downstream boundary conditions for i p becomes 

\p x + iX j \jj = 0 (C-5) 

Similarly, the upstream boundary condition for an outgoing plane wave in the negative 
x-direction is given by 

<^l x -iX^i + M )y>i =0 ( C - 6 ) 

and in terms of \p becomes 

4/ x - iXj^ = 0 (C-7) 

Consider a rectangular region with boundary conditions as shown in figure C-l . 
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\jj = 0 or \}j y + iXj \p = 0 


* = f(y) 



\jj = 0 or + iXji// = 0 


i// = 0 or i//y - iX j i// = 0 
Figure C-1 .—Boundary Conditions for Error Analysis 

We use separation of variables to solve equation (C-3) and obtain the general solution 

2 

i p =^Aj sin px + Bj cos px^ ^2 sin vy+ cos p y ) X| > v 
i // =^A j sinh px + B j ( 
where p and v are related to X] by 


/ w _ 2 2 (C-8) 

1 p =^A] sinh px + B j cosh pxj / A 2 sin vy + B 2 cos U Y ) Xj < v 


2 2 2 
Xj + p + v =0 

2 2 

Consider first the case of X j > v . To satisfy the boundary conditions i// = 0 at y = 0 and 

y = bVK , the parameter v must equal n7r with B-? = 0. 

bVK 


Thus 


where 


imy p -1 

\p = sin • B j cos px + A j sin px 

bVK 


E=-, /X 


2 2 

2 n jt 


b K 


Eigenfunctions for which i// = 0 on all boundaries occur when B j = 0 

m7r 


p = — , m = 0,1, 2,3... 

a 


or 




ntr 


. bVK/ 


(C-9) 


(C-10) 


(C-1 1) 


(C-1 2) 


We take n = 1 and set \p = sin (ny / bVK^ on x = x j . 
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For \p = 0 at x = a, we use 7 ry 

sin — — • sin a - x) 

(C-13) 

sin ju(a ' x i) 

When Xi < — , 
b 


\p = [sin (?ry / bVlC ) • sinh ju( a - x)] / sinh m (a - x j) 


(C-14) 


where 


/i=Vlxi 2 - n V/t 2 K )| 

We now consider the boundary conditions 


(C- 15) 


i// x + i^i 1 / / = 0 atx = a 


(C-16) 


Then equation (C-9) leads to 


7ry 

i// = sin — • • 

b 


fi cos /a(a - x) + iX] sin p( a - x) 


(p cos /i(a - Xj) + iX] sin p(a - xj) j 

We now consider a solution of equation (C-3) with outgoing wave boundary conditions on 
the three boundaries, i.e., 

\p x + iXj \jj = 0 at x = a 


(C-17) 


4*y + iX j i// = 0 at y = 0 
\py -iXj^/ = 0 at y = -Yk b 

In place of the form of the solution given in equation (C-8) we write 


* = 

2 2 2 

where X j + a + q =0 


iay 
c j e + C 2 e 


-loy 


T i 

• dje 


ih(x-a) 


+ d2e 


-iM(x-a)' 


To satisfy the conditions at x = a we obtain for the second bracketed term 

i/r(x-a) -ih(x-a) 

(h - X j ) e +(n + \\)e 


(/a-Xj)e 1M ^ Xl a ^+(M + X])e 
Here, we normalized the factor to make it unity at the upstream boundary x = xj. 


-iM(x r a) 


(C- 18) 


(C-19) 


The constants C| and C 2 are now determined to satisfy the boundary conditions in 
equation (C-18) for the boundaries y = 0 and -b VK . The first boundary condition yields 

(“ + Xj) 

2 ( a - X l) 

while the second boundary condition leads to 
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II 


ia|e 


-iabVK a + X l iabYK" 


a - Ai 


Cl f -iabYK a + x l iabVK" 

- lX l[ e + ¥TXf e J = 0 


or 


. % .2 2 2iabVK 

(a - Xj) - (a + Xj) e =0 

Substituting a = X jz we choose a nonzero root of the equation 

, 2 ^ 2 -2izXibYK 

(z- 1) -Cz+ 1) e 1 = 0 (C-20) 

for each value of Xj. With this specific value of z, the solution to the flat plate equation (C-l) 
is found by multiplying the solution to equation (C-3) by the factor exp(iXjM(x - xj)). 

Thus we obtain 

= e iX j M (x-x j) J JK j zy VK + (. z + l\-iX j zy Vk j 


n 




f 


f(M-Xj) • exp[m(x-a)l +(m + Xj) • exp[- i/u(x - a)] ) ( C -21) 

| (M ~ X i) • exp[iM(xj -a)]+ (ji + Xj) • exp [- i/a (x j - a)] j 


where z is a root of equation (C-20) and 


/a =► iXj VT 


+ z 


In the difference solution, the first bracketed term is used as boundary conditions on the 
mesh line x = X] . Similarly, solutions of the flat plate equation are obtained from equa- 
tions (C-l 3), (C-14), and (C-17). We then have 

V l = |exp[iXjM(x -xj)] • sin(7ry / b) • sinh[/a(a - x)] | / sinh[/u(a - xj)J , Xj < zr / (bVK) 
= jexp[iX 1 M(x-x 1 )] • sin (Try /b) • sin[/i(a-x)] j / sin [/a(a - Xj)] , Xj> 7 r/(bVK) 

(C-22) 


= exp [iXjM (x - x j) J • sin(7ry / b) 


cosh[/a(a - x)] + iXj sinh[/u(a - x)] 

/a cosh j)a(a - Xj)J + iXj sinh[/a(a - Xj)] 


I 


< t r / (bVK) 


= exp[^iXjM(x - Xj)J • sin(7ry / b) • 


j /a cos[/u(a - x)] + iX j sin[/u(a - x)] 
)ju cos jju^a - x j)J + iXj sin[)u(a - 


u)](' 


where /a = 



Xj > 7T / (bVK) 
(C-23) 


The solutions in equations (C-21), (C-22), and (C-23) were used to study the influence of 
mesh boundary conditions on the accuracy of the solutions of the finite difference 
equations. 
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APPENDIX D 


FOURIER ANALYSIS OF TWO STEP FUNCTIONS 
OSCILLATING ON OPPOSITE SIDES OF A PLATE 


Consider a jump in pressure across a plate produced by two pressure step functions 
representing shocks on the two sides of a flat zero thickness airfoil. For the two step 
functions 90° out of phase, the pressure jump is given 

ACp = f(t) = H(x - a sin cot) - H(x - a cos cat) 

where H(x) is the unit function defined by 

H(x) = 0 for x < 0 


(D-l) 


Note that AC p = 0 for x > a and x < -a. This produces a pulse of oscillating width about 


H(x) = 1 for x > 0 
' -a. This produces a 

the range -a < x < a. Furthermore, the center of this pulse also oscillates about x = 0. 

We shall now represent equation (D-l) by a Fourier expansion of the form 

oo oo 

f(t) = ^ b n sin neat + JZ c n c °s neat (D-2) 

n=l n^O 

For 0 < x < a, the coefficients b n are given by 


-1 

•sin (x/a) 

7rb„ = / sin nddd 


-■ 


- TT 

Integration yields 


•- cos (x/a) 
tt - sin (x/a) n 


+ J sin nddd - ^ 


sin n0d0 - f sin nddd (D-3) 
cos' 1 (x/a) 


b n = |(- 1 + cos mr) cos [n sin (x / a)] j / (n7r) 
Replacing sin n d by cos nd in equation (D-3) yields 


(D-4) 


7rc n = j (1 + cos n7r) sin _n sin (x / a) 


-1 


| / n + 2 sin n cos (x / a) J / n (D-5) 


The preceding equations also hold for negative x. Since we are concerned with the ampli- 
tude of the fundamental frequency, we set n = 1 and obtain 


7rb j = - 2 cos |_sin (x / a) 


_1 ' )] = - 2 V 1 - (x/a) 2 


ircj = 2 sin 


[cos (x / a)] = 2l/l - (x / a) 


(D-6) 


(D-7) 
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In the complex form of the amplitude, bj is the imaginary part of exp(io>t), while Cj is the 

real part. Equations (D-6) and (D-7) are plotted in figure 5. Note that the pressure pulse 
is similar to the example in figure 2. 

Thus it appears that the present frequency domain method may well produce the pulse 
amplitude and movement of the fundamental resulting from moving shocks without addi- 
tional treatment such as computing the actual shock motion by a linear analysis and per- 
forming the Fourier analysis. 
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